授業内容(2日目)

• 6.統計解析環境の構築
• 7.データの読み込みと前処理
• 8.記述統計と効果量
• 9.データの可視化
• 10.平均値差と比率差の検定
• 11.回帰分析と因果関係
• 12.因子分析，SEM，メタ分析

SPSSではダメな理由

• 以下の理由からSPSSは，すすめない

1. SPSSでは，できない解析がある

2. 機能のわりに高価（卒業後は使えない）

3. 操作マニュアルは多いが，教育上は悪い影響も（統計学を理解せずに解析するなど）

4. データの整形機能が貧弱

5. 解析過程が残らない（解析が再現出来ない可能性がある）

6. 出力が非常に見にくい（好みの問題ですが・・・）

Excelではダメな理由

• SPSSでデータ整形は難しいので，Excelで前処理をすることも多いのでは？

• しかし，以下の理由から，Excelもすすめない。

1. 操作の過程が残らない（Excelで行った操作を別のファイルに記録しないと再現不可能になる）

2. 生データ自体が改変される

3. 本格的な統計解析・可視化は難しい（HADという選択もあるが・・・）

RとRStuidoを薦める理由

• RとRstudioならSPSSやExcelの問題は，以下のようにクリアできます。

1. できない統計解析はほとんどない

2. 解析過程がしっかり残る

3. データの整形，可視化，解析，解析レポート作成（マニアは論文やプレゼンも）可能

4. こんなにすごいのに，無料！

5. 学生が統計を学習する上で便利（自宅でもできる＆修了後も使える）

Rのインストール

The R Project for Statistical Computingにアクセスし，「download R」をクリックする。

• ミラーサイトを選ぶ(0-CloudかJapanのどれかを選ぶ)

• 使っているOSを選択する

Windowsユーザーの場合

• 「Download R for Windows」をクリックする。

• 「base」をクリックする

• 「Download R-XX for Windows」（xxはバージョン名）をクリックしてダウンロード＆インストールする。

Rtoolsのインストール

• さきほど「base」をクリックしたところで「Rtools」を選択する。

• さきほどインストールしたRのバージョンと同じバージョンのRtoolsをクリックする。

• 「RtoolsXX installer」(XXはバージョン)をクリックして，ダウンロード＆インストールする。

Windowsユーザーの注意事項

ユーザー名を日本語名にしていると問題が発生することがある。英語名のユーザーを追加するか，こちらを参考にして，設定を変更する。

Macユーザーの場合

Xcodeのインストール

• まず，「Xcode」をインストールする。App Storeを起動して，「Xcode」で検索して，「入手」する（「開く」になっている場合は，インストール済みである）。インストールには時間かかる。

コマンドラインツールのインストール

• インストール後に以下をターミナルで実行する。

xcode-select --install

Xquartzのインストール

• Xquartzにアクセスし，ダウンロード＆インストールする。

Rのインストール

• 「Download R for mac OS」をクリックする。

• Macの種類に合わせてクリックして，ダウンロード＆インストールする。Apple siliconの場合はarm64版をインストールし，Intel Macの場合はx86_64版をインストールする。

RStuidoのインストール

• posit社のサイトにアクセスして，「DOWNLOAD RSTUDIO」をクリックする。

• 移動した先のサイトで下にスクロールして，Free版の「RStudio Desktop」の「DOWNLOAD」をクリックする。

• 移動した先のサイトで下にスクロールして，自分の使っているOSのものをダウンロードしてインストールする（Macはアプリケーションにコピーする）。

Rの関数：標準関数

• Rは標準関数で色々な統計解析ができる。

• Positチートシートの「Base R」（日本語版）を参照したり，「R　使いたい解析名」などでグーグル検索して調べる。

• Rの関数名はわかっているけど，どう使ったらわからない場合，「?関数名」をコンソールにタイプするとヘルプがでてくる。

Rと関数

• Rを使いこなすには，関数の理解が大切です。以下のように，関数の（）内で，引数を指定すると，左辺に結果が代入されます。

戻り値（結果） = 関数（引数（入力））

例) 平均値の関数mean()の場合

結果 = mean(平均したいデータ)

Rの使い方：CRANパッケージ

• Rでは，標準関数以外にも様々な解析パッケージが作成され，多くはCRANに登録されている。

• CRANにあるパッケージは以下のようにコンソールにタイプするとインストールできます。

install.packages("パッケージ名", dependencies = TRUE)

• パッケージを使う場合は，以下のようにコンソールにタイプして，呼び出します。

library(パッケージ名)

パッケージのインストール

今回の講義で使用するパッケージをインストールする。

install.packages(c('tidyverse','netmeta','meta',
                   'esc','igraph','lme4','lmerTest',
                   'phia','devtools','lavaan',
                   'gtsummary','effsize','ggdist',
                   'corrr','coin','pwr','psych','paran'))

dmetarはgithub経由でインストールする。

devtools::install_github("MathiasHarrer/dmetar")

パッケージや関数の使い方が分からない

• パッケージを探す：「R 解析名」でグーグル検索

• 関数の使い方を探す：（１）「?関数名」をRStudioのコンソールにタイプするとヘルプがでてくるので読む，（２）CRANでパッケージを検索してマニュアルを読む（分かりやすいビネットが公開されていることもある），（３）「関数名　解析名」などでグーグル検索（英語のほうが情報が多そうだがRは日本語でも重要な情報が見つかることも多い）

パッケージや関数の使い方が分からない

• プログラミングならchatGPTが便利です。

Rはエラーが出て怖い！

• エラー内容（赤字）をグーグル検索する or chatGPTに聞いてみる

Quarto

• Consoleにコマンドを打ち込むやり方だと手元に解析過程が残らない，明示的にコマンドとその説明文章を残す習慣をつける必要がある。

• Quartoを使うと，便利なMarkdown記法を使って，文章と解析コードを残せる。

• Quartoは複数の出力(HTML, PDF, Word)ができ，R以外の言語も対応し，RStudio以外でも利用できる（Quartoは元々あったRmarkdownを統合したものという位置づけ）。

Quarto Documentの作成

• 新規作成から，Quarto Documentを選択する（ちなみに，この講義資料はQuarto Presentationを使用している）。

• Titleを適当につけてCreateをクリックする。

Quartoに書き込む

• SourceモードとVisualモードがあるが，VisualモードはWordみたいな感じで作業ができる。

• ①設定情報。②本文（見出しや引用や画像の挿入）。③解析部分のチャンク。④RenderをクリックするとHTML出力される。

ワーク

• RとRStudioとRとパッケージをインストールして，Quartoドキュメントを作成して（書き込みはしなくてもOKです），Renderしてみよう！

• 時間が空いたなら，以下の資料を読んで再現可能な環境構築について学ぼう。

再現可能な日本語論文執筆入門:jpaRmdで実現する再現可能で低コストな日本語論文執筆のはじめの一歩

統計学とプログラミングの学習

• 統計学やプログラミングは実際に手を動かして学ぶのが効率的

• Rのパッケージには統計の学習用のデータも用意されていますが，キレイ整いすぎている。そこで，実際のオープンデータを使うのをオススメする。オープンデータを使った解析は，２次分析研究にもなり得る。

• ２次分析(Secondary Data Analysis)とは，既存データ（preexisting data）に，既存の研究とは違う目的（視点）をもって分析することである。

既存データの種類(Weston et al., 2019)

1. 大規模調査のデータ（政府統計，パネル調査）

2. 個別の研究データ（各研究室のデータ）

3. ビッグデータ（ウェブスクレイピング，医療報酬の明細書など）

• 既存データの中にはオープンデータも増えている。

• オープンデータとは，「自由に使えるデータ」（庄司, 2014)

→どのようなライセンスの下でデータが公開されているのかが重要になる

オープンデータはどこに？

オープンデータのリストとデータ検索サイト

• 心理学関連オープンデータリスト（国里と有志の方がまとめたもの）

• Open Science Framework

• Google Dataset Search

• DataCite Search

• Data Citation Index (Web of Scienceの一部,有料）

Data Journal

データを掲載するジャーナルも増えてきている。

• Journal of Open Psychology Data (JOPD)

• Scientific Data

• Data in Brief

国里も参加した新型コロナ禍中の国際共同プロジェクトのデータはScientific Dataに掲載（世界で73,223名が参加した研究）

Buchanan, E.M., Lewis, S.C., Paris, B. et al. The Psychological Science Accelerator’s COVID-19 rapid-response dataset. Sci Data 10, 87 (2023).

今回使用するデータ

• 以下の論文のデータを使う。

Woodworth, R. J., O’Brien-Malone, A., Diamond, M. R., & Schüz, B. (2018). Data from, “Web-based Positive Psychology Interventions: A Reexamination of Effectiveness.” Journal of Open Psychology Data. https://doi.org/10.5334/jopd.35

• データは以下のfigshareで公開されている。

https://figshare.com/articles/dataset/A_randomized_placebo-controlled_trial_of_positive_psychology_interventions_in_Australia/1577563/1

Woodworth et al.(2018)のデータ

対象者：一般参加者295名

介入：３つのポジティブ心理学的介入と統制条件に無作為割付。(1)“Using Signature Strengths（ウェブサイトで性格の強みを評価し，上位の強みのうちの１つを次の週に新しい別のやり方で実施する）”，(2)“Three Good Things(その日に起こった良いことを３つあげて，なぜそれが起こったのか原因説明とともに記録する)”, (3)“Gratitude Visit(親切してくれたけど感謝を示してない人に手紙を書く介入)”,(4)“Recording early memories(統制条件，幼いことを記憶を書く)”

アウトカム：Authentic Happiness Inventory (AHI) の24項目，Center for Epidemiological Studies Depression(CES-D) scaleの20項目

測定時期：介入前(0), 介入後(1), 1週間後(2), １ヶ月後(3)，３ヶ月後(4)，６ヶ月後(5)

Woodworth et al.(2018)のデータにアクセスして，ファイルをダウンロードして開きます。

RStudioでプロジェクトを作成

• RStudioにはプロジェクトという機能があり，フォルダやファイルやパッケージの管理に便利です。研究プロジェクト単位でプロジェクトを作って，関連するファイルをいれると解析に便利です。

• FileからNew Projectをクリック

• New Directoryをクリック

• New Projectをクリック

• Directory名を英数字＋アンダーバーでつける。Browseをクリックしてディレクトリを作る場所を指定する（今はいいけど，将来的には”Create a git repository”や”Use renv with this project”も使いこなしてほしい）。

• 作成したwoodworth_2018というディレクトリ内に，“data”という名前のディレクトリを作成して，先程ダウンロードしたWoodworth et al.(2018)のデータのファイルをいれてください。

• これで下準備完了です。以降はFile > Recent Projectsをみると最近使ったプロジェクトがでてくるので，それをクリックすれば，すぐに作業環境に戻れます（Recent Projectsに出てこない場合は，Open Projectをクリックして，当該ディレクトリに移動して，XXX.Projというファイルを開いてください）。

Quartoドキュメントの準備

• Quartoドキュメントを用意します。新規作成ボタンからQuarto Documentをクリックして，タイトルを入れます。なお，タイトルは何でもいいです。

• ファイルを開いたらSaveボタンをクリックして，ファイルに名前をつけて保存します。その際に，ファイル名は，analysis_01, analysis_02と順番が分かる名前にしておくと良いです（今回はanalysis_01にしましょう）。

• こんな感じになります。なお，最初に開いたファイルは日本語入力できないという謎バグがあるので，適当なファイルを開いて保存してから（私はText fileを新規作成して，README.mdって名前で保存します），再度analysis_01.qmdを開きます。

• 設定部分より下は削除して，日本語入力してみましょう。書式設定とかできるので，適当にいじってみましょう。また，見出しはとても重要です。文章の構造化を意識しましょう。

• Rチャンクの追加は，Insert -> Code Cell -> Rでできます。

• Rチャンクが用意できました。以降の作業は，このanalysis_01.qmdに書き込んでいきましょう。

データの読み込み

• データの読み込み，データの前処理，可視化などはtidyverseパッケージ使うと便利です。

• dataフォルダに保存したWoodworth et al.(2018)のデータを読み込みます。データはcsv形式だったので，read_csv()で読み込みます。ファイルはparticipant-info.csvとahi-cesd.csvの２種類があります。dataフォルダにいれているので，ファイル名の前に”data/“をいれています。こうすると，dataフォルダ内の特定のファイルを指定できます。

library("tidyverse") 
participant <- read_csv("data/participant-info.csv")
data <- read_csv("data/ahi-cesd.csv") 

• qmdだと以下のような感じです。Rチャンク右端の緑の三角形をクリックするとコードが実行されますので，クリックしてください。

• コードを実行するとEnvironmentにオブジェクトが保存されます。右端のアイコンをクリックするとオブジェクト（今回はデータ）を見ることができます。

• participant-info.csvには，id, intervention, sex, age, educ, incomeのデータが入っています。

• ahi-cesd.csvには，id, occasion, elapsed.days, intervention, ah01…などのデータが入っています。

• 個人的にはエクセルファイルはcsv形式に変換して処理するのが良いですが，エクセルから直接データを読みたい場合は，readxlパッケージのread_excel()を使って読み込むことができます。
• foreignパッケージを使うと，‘Epi Info’, ‘Minitab’, ‘S’, ‘SAS’, ‘SPSS’, ‘Stata’, ‘Systat’, ’Weka’などのファイルも読み込めます。

ネイティブパイプの使用

• Rではパイプ演算子という便利なものを使います。以前はmagrittrパッケージを用いた%>%だったのですが，R4.1からは標準で|> が使えます。

• |> が使えるように設定します。Tools -> Global options -> Codeで，Use native pipe operator…にチェックをいれてApplyする。

Windowsの場合Ctrl + Shift + m　Macの場合Command + Shift + mのショートカットでパイプ演算子が入力できます。

データフレーム

• Rではデータフレームという構造でデータを扱う。データフレームは，a行b列からなる表形式のデータであり，列には変数が含まれる。

• read_csv()で読み込んだデータはデータフレームになっており，以下のオブジェクト名を打つとデータフレームが表示される。

• Rでの代表的なデータ型としては以下がある（read_csv()で読み込むと自動で判定をする）。

データの前処理

• 前処理にはtidyverseパッケージに含まれているdplyrパッケージを使う。色々と機能はあるが，ミニマムには以下の８つを覚える。

• filter() 行の絞り込み

• arrange() 行の並び替え

• select() 指定した列の抽出

• mutate() 列の追加

• group_by() グループ化

• summarize() 要約

• pivot_longer() 横データから縦データへ変換（時系列解析などは縦データが便利）

• pivot_wider() 縦データから横データへ変換

filter()

• filter() で行の絞り込みができます。性別を女性(1,男性は2)だけに絞り込む場合は以下のようにできます。パイプ演算子は，左にあるオブジェクトを右（改行していることが多い）の関数にいれる機能を持っています。パイプ演算子を使うと処理を重ねていけるので見やすいし便利です。今回は，participantデータをfilter()にいれて，sexが1のものだけを抽出しています。

• パイプ演算子の挿入はWindowsの場合Ctrl + Shift + m　Macの場合Command + Shift + mです。

participant |> 
  filter(sex == 1)

# A tibble: 251 × 6
      id intervention   sex   age  educ income
   <dbl>        <dbl> <dbl> <dbl> <dbl>  <dbl>
 1     2            1     1    59     1      1
 2     3            4     1    51     4      3
 3     4            3     1    50     5      2
 4     6            1     1    31     5      1
 5     7            3     1    44     5      2
 6     8            2     1    57     4      2
 7     9            1     1    36     4      3
 8    10            2     1    45     4      3
 9    11            2     1    56     5      1
10    12            2     1    46     4      3
# ℹ 241 more rows

• filter()では，XX以上でしぼりこみもできる。以下は，年齢が40以上で絞り込んでいる。summarize() を使うと，変数名=関数（変数）で要約ができます。最小40歳，最大85歳ですね。

participant |> 
  filter(age >= 40) |> 
  summarize(min_age = min(age), max_age = max(age))

# A tibble: 1 × 2
  min_age max_age
    <dbl>   <dbl>
1      40      83

以下のように４０歳以下かつ女性で絞り込むこともできます。

participant |> 
  filter(age <= 40 & sex == 2)

# A tibble: 17 × 6
      id intervention   sex   age  educ income
   <dbl>        <dbl> <dbl> <dbl> <dbl>  <dbl>
 1     1            4     2    35     5      3
 2    15            4     2    27     2      2
 3    44            2     2    36     4      3
 4   130            4     2    21     1      1
 5   140            3     2    31     3      1
 6   193            2     2    25     5      2
 7   204            4     2    32     5      1
 8   210            1     2    28     4      2
 9   234            3     2    33     5      3
10   235            2     2    18     2      1
11   237            4     2    33     3      3
12   243            4     2    26     4      2
13   251            4     2    31     5      2
14   263            1     2    32     5      2
15   269            3     2    23     4      1
16   284            3     2    19     4      1
17   288            3     2    26     5      1

arrange()

• arrange() で行の並び替えができます。指定した変数で昇順で並び替えます。

participant |> 
  arrange(age)

# A tibble: 295 × 6
      id intervention   sex   age  educ income
   <dbl>        <dbl> <dbl> <dbl> <dbl>  <dbl>
 1   235            2     2    18     2      1
 2   267            4     1    19     2      1
 3   284            3     2    19     4      1
 4    20            3     1    21     4      1
 5   130            4     2    21     1      1
 6   230            2     1    21     2      1
 7    23            1     1    22     4      1
 8   122            1     1    23     5      2
 9   138            1     1    23     2      1
10   247            1     1    23     4      1
# ℹ 285 more rows

• arrange(-age)のように-を使うことで降順にすることもできます。

participant |> 
  arrange(-age)

# A tibble: 295 × 6
      id intervention   sex   age  educ income
   <dbl>        <dbl> <dbl> <dbl> <dbl>  <dbl>
 1    51            2     1    83     2      2
 2   244            2     1    77     3      2
 3   215            4     1    75     3      2
 4    83            4     1    71     2      1
 5   114            1     1    71     4      1
 6   155            3     2    71     4      1
 7    93            3     1    68     1      2
 8   141            2     2    67     5      2
 9   281            4     2    67     4      2
10    67            1     1    65     4      2
# ℹ 285 more rows

select()

• select()で指定した列の抽出ができます。例えば，id,age,sexだけ抽出したい場合は，以下のように選択できます。

participant |> 
  select(id,age,sex)

# A tibble: 295 × 3
      id   age   sex
   <dbl> <dbl> <dbl>
 1     1    35     2
 2     2    59     1
 3     3    51     1
 4     4    50     1
 5     5    58     2
 6     6    31     1
 7     7    44     1
 8     8    57     1
 9     9    36     1
10    10    45     1
# ℹ 285 more rows

mutate()

• mutate() で列の追加ができます。例えば，dataにはCES-Dの得点が入っていますが，その中でポジティブ感情の項目（4,8,12,16）だけを合計して得点化します。確認のためにsummarizeもしています。

data |> 
  mutate(positive_score = cesd04 + cesd08 + cesd12 + cesd16) |> 
  summarize(mean = mean(positive_score),
            sd = sd(positive_score),
            min = min(positive_score),
            max = max(positive_score))

# A tibble: 1 × 4
   mean    sd   min   max
  <dbl> <dbl> <dbl> <dbl>
1  12.2  3.28     4    16

group_by()

• group_by()でグループ化します。dataのCES-Dの合計得点(cesdTotal)の介入前の得点(occasionが0)に対して，interventionでグループ分けして，summarizeします。

data |> 
  filter(occasion == 0) |> 
  select(intervention, cesdTotal) |> 
  group_by(intervention) |>
  summarize(mean = mean(cesdTotal),
            sd = sd(cesdTotal),
            n = n())

# A tibble: 4 × 4
  intervention  mean    sd     n
         <dbl> <dbl> <dbl> <int>
1            1  15.1 11.7     72
2            2  16.2  9.89    76
3            3  16.1 11.8     74
4            4  12.8  9.51    73

pivot_longer()とpivot_wider()

• pivot_longer()で横データから縦データへ変換，pivot_wider()で縦データから横データへ変換します。dataはすでにlong形式になっています（occasionごとに横にデータが並んでいるのではなく，１つの参加者でoccationごとに複数行が並んでいます）。これを横データにしてみます。select(id,occasion,intervention,cesdTotal)で少し変数をしぼって，pivot_wider(names_from = occasion, values_from = cesdTotal) で横データにします。最後に，-> data_wideとするとdata_wideというオブジェクト名として保存されます。occasionの１以降は回答してない人も多いので，NULLになっていますね。

• idの8がoccasionの2で二回回答しているので，二回目の回答を除外します(指定がしにくくてid とcesdTotalで指定しています)。

data |> 
  filter(!(id == 8 & cesdTotal == 49)) |> 
  select(id,occasion,intervention,cesdTotal) |> 
  pivot_wider(names_from = occasion, values_from = cesdTotal) -> data_wide

• 今度は，pivot_longer()で縦データを横データに変換します。pivot_longer(cols = c("移動する列"), names_to = "分ける変数名", values_to = "値の名前") で横データになります。ただ，元に戻るのではなく，横→縦変換時に欠測になっているので，横データに戻した時に欠測が含まれています。

data_wide |> 
  pivot_longer(cols = c(-id,-intervention), names_to = "occasion", values_to = "cesdTotal")

# A tibble: 1,770 × 4
      id intervention occasion cesdTotal
   <dbl>        <dbl> <chr>    <list>   
 1     1            4 0        <dbl [1]>
 2     1            4 1        <dbl [1]>
 3     1            4 2        <NULL>   
 4     1            4 3        <NULL>   
 5     1            4 4        <NULL>   
 6     1            4 5        <NULL>   
 7     2            1 0        <dbl [1]>
 8     2            1 1        <dbl [1]>
 9     2            1 2        <dbl [1]>
10     2            1 3        <dbl [1]>
# ℹ 1,760 more rows

ワーク

Woodworth et al.(2018)のデータについて，８つの関数を駆使して，様々にデータを変形してみましょう！

glimpse()

• Environmentでも確認できますが，glimpse()でデータの変数を確認できる。

participant |> 
  glimpse()

Rows: 295
Columns: 6
$ id           <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
$ intervention <dbl> 4, 1, 4, 3, 2, 1, 3, 2, 1, 2, 2, 2, 4, 4, 4, 4, 3, 2, 1, …
$ sex          <dbl> 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, …
$ age          <dbl> 35, 59, 51, 50, 58, 31, 44, 57, 36, 45, 56, 46, 34, 41, 2…
$ educ         <dbl> 5, 1, 4, 5, 5, 5, 5, 4, 4, 4, 5, 4, 5, 1, 2, 1, 4, 5, 3, …
$ income       <dbl> 3, 1, 3, 2, 2, 1, 2, 2, 3, 3, 1, 3, 3, 2, 2, 1, 2, 2, 1, …

summary()

• 連続変数に関しては，summary()で必要な記述統計を出すことができる（selectを組み合わせると便利）。

participant |> 
  select(age, educ, income) |> 
  summary()

      age             educ          income     
 Min.   :18.00   Min.   :1.00   Min.   :1.000  
 1st Qu.:34.00   1st Qu.:3.50   1st Qu.:2.000  
 Median :44.00   Median :4.00   Median :2.000  
 Mean   :43.76   Mean   :3.98   Mean   :2.044  
 3rd Qu.:52.00   3rd Qu.:5.00   3rd Qu.:3.000  
 Max.   :83.00   Max.   :5.00   Max.   :3.000  

table()

• 質的変数の場合は，table()でクロス集計を出すことができる。このデータは女性が多いことが分かる。

participant |> 
  select(sex) |> 
  table()

sex
  1   2 
251  44 

gtsummary()その１

• gtsummary パッケージのtbl_summary()を使うと簡単に変数の要約ができます。

library(gtsummary)
participant |>
  select(-id) |> 
  tbl_summary()

gtsummary()その２

• byでグループ分けをして表を作ることもできる。

participant |>
  select(-id) |> 
  tbl_summary(by = intervention)

• gtsummary()の調整法としては，こちらが参考になります。

• より柔軟に表を作成するパッケージとしてはgtパッケージがありますが，説明を省略します。

連続変数における２群の差の効果量

２群の差は，標準化平均値差(Hedgesのg)を用いる

• \(Hedgesのg = \frac{介入群の平均-統制群の平均}{２群をプールした標準偏差}\)

• \(2群をプールした標準偏差=\sqrt{\frac{(n_{介入}-1)\hat{\sigma}^{2}_{介入}+(n_{統制}-1)\hat{\sigma}^{2}_{統制}}{(n_{介入}+n_{統制}-2)}}\)

\(\hat{\sigma}^{2}\)は不偏分散。他にCohenのdやGlassの⊿もあるが，バイアスを修正した効果量のHedgesのgがよく使用される

• \(標準化平均値差の分散=\frac{n_{介入}+n_{統制}}{n_{介入}n_{統制}}+\frac{g^2}{2(n_{介入}+n_{統制})}\)

• \(標準化平均値差の標準誤差=\sqrt{標準化平均値差の分散}\)

分散と標準偏差は効果量とサンプルサイズで決まる

標準化平均値差(Hedgesのg)を計算する

• dataの介入前の(3)“Gratitude Visit(親切してくれたけど感謝を示してない人に手紙を書く介入)”と(4)“Recording early memories(統制条件，幼いことを記憶を書く)”を比較してみる。3の方が得点が高い。なお，表はgtパッケージのgt()できれいにした。

library(effsize)
library(gt)

data |> 
  filter(intervention == 3 | intervention == 4) |> 
  filter(occasion == 0) |> 
  select(intervention,cesdTotal) |> 
  group_by(intervention) |> 
  summarize(mean = mean(cesdTotal),sd = sd(cesdTotal), n = n()) |> 
  gt()

Hedgesのgを計算すると弱い効果量がある。

Cohen's d

d estimate: 0.3010943 (small)
95 percent confidence interval:
     lower      upper 
-0.0267866  0.6289752 

２つの連続変数の関連の効果量

２変数間の関連の強さは，ピアソンの積率相関係数を用いる

• \(r = \frac{共分散_{x,y}}{標準偏差_{x}標準偏差_{y}}\)

• \(相関係数の分散 = \frac{(1-r^2)^2}{n-1}\),\(相関係数のSE = \frac{1-r^2}{\sqrt{n-1}}\)

• 統計学的特性の望ましさから(サンプルサイズが小さいとrのseにはバイアスが生じる)，フィッシャーのZ変換を用いることが多い

• \(z = 0.5ln(\frac{1+r}{1-r})\), \(zの分散 = \frac{1}{n-3}\), \(zのSE = \frac{1}{\sqrt{n-3}}\)

• zからrへは\(\frac{exp(2z)-1}{exp(2z)+1}\)で変換できる(信頼区間も同様)。

２値変数（あり・なし）の効果量

1. リスク差=介入群のありの比率-統制群のありの比率

2. \(治療必要数=\frac{1}{リスク差}\)

3. \(リスク比=\frac{介入群のありの比率}{統制群のありの比率}\)

4. \(オッズ比=\frac{\frac{介入群のありの比率}{1-介入群のありの比率}}{\frac{統制群のありの比率}{1-統制群のありの比率}}\)

リスク差や治療必要数(Number needed to treat)は解釈しやすいが，計算上オッズ比が便利なので，オッズが使用される。さらに対数オッズ比(ln(オッズ比))を使う（1を起点に解釈→0を起点に解釈，分散も計算がしやすい）。

オッズ比の分散=介入群あり比率+介入群なし比率+統制群あり比率+統制群なし比率

ワーク

Woodworth et al.(2018)のデータについて，論文のTable1に該当するような表を作成してみよう！他の群間の平均値を比べて，効果量を計算してみよう！

データは必ず可視化しよう！

• データに対してすぐにでも高度な統計手法を適用したくなるが，まずは記述統計とともにデータの可視化を必ずしましょう。
• データの可視化は，(1)データに対する前処理が適切にできているかチェック，(2)データの測定法が適切にできているかチェック，(3)データの傾向をチェック，(4)データに適用する統計手法が適切かチェックするために行う。
• データの可視化を行ってないとエラーに気づかない可能性があり，結果を歪めてしまう可能性がある。

ggplot2

• Rには可視化のための関数があるが，ggplot2パッケージがきれいに図が作成できて，コードも読みやすく書ける。

• ggplot2は，グラフィックの文法 (grammar of graphics)という思想の元に作成されている。これは，有る決まった構造で作成されたレイヤーを重ねていって図を作成するという考え方です。

棒グラフ

• 棒グラフを題材にggplot2の使い方を学びましょう。ggplot2はtidyverseパッケージを読み込むと使えます。まずは，性別ごとの人数の棒グラフを描きます。女性か男性か図に書き込みたいので，sexの1,2の数値を女性と男性に変換した変数をsex2として追加します。

library(tidyverse)
participant |> 
  mutate(sex2 = ifelse(sex ==1, '女性','男性')) -> participant

• まずggplotに使用するデータを指定します。この段階では図を書くキャンバスが出力されるだけです。

ggplot(data = participant) 

• キャンバスに追加します。追加する際にはパイプ演算子ではなくて+を使います。aes()では，x軸, y軸の変数や色などを指定します。ここでは，１変数のsex2だけ指定します。macの場合は文字化けしているのですが，ちょっと置いておきます。

ggplot(data = participant) +
  aes(x=sex2)

• 続けて，geom_bar()で，棒グラフを重ねます。

ggplot(data = participant) +
  aes(x=sex2) +
  geom_bar()

• 文字化けが気になってきたので，plotのテーマはtheme_grayのままbase_family = "HiraKakuPro-W3" を指定します。

ggplot(data = participant) +
  aes(x=sex2) +
  geom_bar() +
  theme_gray(base_family = "HiraKakuPro-W3")

• ちょっと灰色の背景だと論文に載せにくいなと思ったかもしれません。その場合は，theme_classic()を使います。あれ？グラフが浮いている感じして気持ち悪いですね。

ggplot(data = participant) +
  aes(x=sex2) +
  geom_bar() +
  theme_classic(base_family = "HiraKakuPro-W3") 

• scale_y_continuous(expand=c(0,0)) を追加すると浮かなくなります。

ggplot(data = participant) +
  aes(x=sex2) +
  geom_bar() +
  theme_classic(base_family = "HiraKakuPro-W3") +
  scale_y_continuous(expand=c(0,0))

• せっかくなので，男女でグラフの色を分けたいと思います。aes()内のfillにsex2を指定します。

ggplot(data = participant) +
  aes(x=sex2,fill=sex2) +
  geom_bar() +
  theme_classic(base_family = "HiraKakuPro-W3") +
  scale_y_continuous(expand=c(0,0))

• ちょっと色が派手で灰色と白くらいにしたいです。

ggplot(data = participant) +
  aes(x=sex2,fill=sex2) +
  geom_bar() +
  theme_classic(base_family = "HiraKakuPro-W3") +
  scale_y_continuous(expand=c(0,0)) +
  scale_fill_manual(values=c("white","gray"))

あ！線がない！geom_bar()内に色を追加します。

ggplot(data = participant) +
  aes(x=sex2,fill=sex2) +
  geom_bar(color="black") +
  theme_classic(base_family = "HiraKakuPro-W3") +
  scale_y_continuous(expand=c(0,0)) +
  scale_fill_manual(values=c("white","gray"))

• 図はいいですが，x軸とy軸のラベルを変えたいですね。labs()を追加します。

ggplot(data = participant) +
  aes(x=sex2,fill=sex2) +
  geom_bar(color="black") +
  theme_classic(base_family = "HiraKakuPro-W3") +
  scale_y_continuous(expand=c(0,0)) +
  scale_fill_manual(values=c("white","gray")) +
  labs(x = "性別", y = "人数")

• この図の場合は，凡例は不要ですね。theme(legend.position = "none")を追加します。

ggplot(data = participant) +
  aes(x=sex2,fill=sex2) +
  geom_bar(color="black") +
  theme_classic(base_family = "HiraKakuPro-W3") +
  scale_y_continuous(expand=c(0,0)) +
  scale_fill_manual(values=c("white","gray")) +
  labs(x = "性別", y = "人数") +
  theme(legend.position = "none")

使えそうな図が作れました。ggplot2は面倒な感じもしなくもないですが，キャンバスにちょっとずつレイヤーを足していくことで図を作る感じをイメージするといいです。その都度，「ggplot 凡例　消す」とかで検索すると簡単に情報が手に入ると思います（chatGPTも便利です）。

ヒストグラム

• 量的な変数についてはヒストグラムを書いて分布を確認します。

• dataからoccasion=0のCES-D得点のヒストグラムを描きます。filterでoccasion=0に絞り込んで，それをggplotに投げます。パイプ演算子を使っているので，filterされたdataがggplotに入るので，aesだけ指定します。geom_histogram()でヒストグラムが描けます。

data |> 
  filter(occasion == 0) |> 
  ggplot(aes(x=cesdTotal)) +
  geom_histogram()

• theme_classic()にしたり，グラフが浮かないようにしたり，x軸のラベルを変えたりして，調整します。

data |> 
  filter(occasion == 0) |> 
  ggplot(aes(x=cesdTotal)) +
  geom_histogram() +
  theme_classic() +
  scale_y_continuous(expand=c(0,0)) +
  labs(x = "CES-D") 

• ヒストグラムも良いですが，密度プロット（density plot）もあると便利です。ヒストグラムと重ねて描くために，geom_histogram(aes(y=after_stat(density)),fill='white', color='black') という指定をしています。また，geom_density(fill='black', alpha = 0.5) のalphaは重ねた際の透過度になります。

data |> 
  filter(occasion == 0) |> 
  ggplot(aes(x=cesdTotal)) +
  geom_histogram(aes(y=after_stat(density)),fill='white', color='black') +
  geom_density(fill='black', alpha = 0.5) +
  theme_classic() +
  scale_y_continuous(expand=c(0,0)) +
  labs(x = "CES-D", y="") 

散布図

• dataからahiTotalとcesdTotalの散布図を描いてみます。相関や回帰の前には散布図を描くようにしましょう。散布図はgeom_point()で描けます。

data |> 
  filter(occasion == 0) |> 
  ggplot(aes(x = cesdTotal, y = ahiTotal)) +
  geom_point()

• 論文用に調整します。きれいに負の相関をしています。

data |> 
  filter(occasion == 0) |> 
  ggplot(aes(x = cesdTotal, y = ahiTotal)) +
  geom_point() +
  theme_classic() +
  labs(x = "CES-D", y = "AHI")

• 散布図内にサブグループがある場合に，色を塗り分けたいこともあるかもしれない。fillとcolorに群の情報をいれると塗り分けができる。なお，interventionは数値型になっているので因子型に変換している。

data |> 
  filter(occasion == 0) |> 
  mutate(intervention = as.factor(intervention)) |> 
  ggplot(aes(x = cesdTotal, y = ahiTotal, 
             fill = intervention, color = intervention)) +
  geom_point() +
  theme_classic() +
  labs(x = "CES-D", y = "AHI")

• ３つの変数の関係をみる上では，バブルプロットも便利です。sizeに３つ目の変数を指定すれば描けます。今回はあまり意味のあるものではないので，そこから意味は読み取れないですね。

data |> 
  filter(occasion == 0) |> 
  mutate(intervention = as.factor(intervention)) |> 
  ggplot(aes(x = cesdTotal, y = ahiTotal, size = intervention)) +
  geom_point() +
  theme_classic() +
  labs(x = "CES-D", y = "AHI")

折れ線

• 時間的な変化や回帰直線などをプロットする際には，折れ線が使える。

• 今回は，occasionが0~5まで５時点あるので，id=1〜10までプロットしてみる。工夫としては，idで10以下のデータにしぼったこと，idは因子型にしたことになる。折れ線は，geom_line()を使う。

data |> 
  filter(id <= 10) |> 
  mutate(id = as.factor(id)) |> 
  ggplot(aes(x = occasion, y = cesdTotal, color = id)) +
  geom_line() 

• theme_classic()などの見栄えの調整をした。geom_point()も重ねると点がわかりやすいので追加している。

data |> 
  filter(id <= 10) |> 
  mutate(id = as.factor(id)) |> 
  ggplot(aes(x = occasion, y = cesdTotal, color = id)) +
  geom_line() +
  geom_point() +
  theme_classic() +
  labs(x = "Time", y = "CES-D")

箱ひげ図

• 群間の比較では箱ひげ図が便利と思われる。箱ひげ図は，最小，第１四分位(下から25%)，中央値，第３四分位(下から75%)，最大値，外れ値をプロットする。
• dataの介入前の(3)“Gratitude Visit(親切してくれたけど感謝を示してない人に手紙を書く介入)”と(4)“Recording early memories(統制条件，幼いことを記憶を書く)”を比較してみる。箱ひげ図はgeom_boxplot()でプロットできる（stat_boxplot()はエラーバーを出すため）。

data |> 
  filter(intervention == 3 | intervention == 4) |> 
  filter(occasion == 0) |> 
  mutate(intervention = as.factor(intervention)) |> 
  ggplot(aes(x = intervention, y = cesdTotal)) +
  stat_boxplot(geom = "errorbar", width = 0.3) +
  geom_boxplot()

• 介入名を入れたり，見栄えの調整をする。

data |> 
  filter(intervention == 3 | intervention == 4) |> 
  filter(occasion == 0) |> 
  mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |> 
  ggplot(aes(x = intervention2, y = cesdTotal)) +
  stat_boxplot(geom = "errorbar", width = 0.3) +
  geom_boxplot() +
  theme_classic() +
  labs(x = "Intervention", y = "CES-D")

• 箱ひげ図よりもデータの分布をみたい場合は，バイオリンプロットが便利です。geom_violin(trim = FALSE)で描けます。

data |> 
  filter(intervention == 3 | intervention == 4) |> 
  filter(occasion == 0) |> 
  mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |> 
  ggplot(aes(x = intervention2, y = cesdTotal)) +
  geom_violin(trim = FALSE) +
  theme_classic() +
  labs(x = "Intervention", y = "CES-D")

• バイオリンプロットにさらにデータポイントを描く（ドットプロット）とよりイメージしやすくなります。

data |> 
  filter(intervention == 3 | intervention == 4) |> 
  filter(occasion == 0) |> 
  mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |> 
  ggplot(aes(x = intervention2, y = cesdTotal)) +
  geom_violin(trim = FALSE) +
  geom_jitter(position=position_jitter(0.1))+
  theme_classic() +
  labs(x = "Intervention", y = "CES-D")

レインクラウドプロット(Raincloud plot)

• 箱ひげ図，バイオリンプロット，ドットプロットを組み合わせたものとして，レインクラウドプロットがある。ggdistパッケージで描けます。まず，stat_halfeye()でバイオリンプロットの半分みたいなものを描きます。

library(ggdist)
data |> 
  filter(intervention == 3 | intervention == 4) |> 
  filter(occasion == 0) |> 
  mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |> 
  ggplot(aes(x = intervention2, y = cesdTotal)) +
  stat_halfeye()

• ドットプロット追加します。stat_halfeye()の信頼区間っぽいやつを削除するために，point_color = NA, .width = 0を指定します。

library(ggdist)
data |> 
  filter(intervention == 3 | intervention == 4) |> 
  filter(occasion == 0) |> 
  mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |> 
  ggplot(aes(x = intervention2, y = cesdTotal)) +
  stat_halfeye(width = .5, point_color = NA, .width = 0) +
  stat_dots(side = "left", dotsize = 1, justification = 1.05, binwidth = 1)

• 箱ひげ図を追加します。

library(ggdist)
data |> 
  filter(intervention == 3 | intervention == 4) |> 
  filter(occasion == 0) |> 
  mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |> 
  ggplot(aes(x = intervention2, y = cesdTotal)) +
  stat_halfeye(width = .5, justification = -.15,point_color = NA, .width = 0) +
  stat_dots(side = "left", dotsize = 1, justification = 1.1, binwidth = 1) +
  geom_boxplot(width = .1,outlier.shape = NA) 

• 最終調整をします。なぜか左に寄っちゃうので，coord_cartesian(clip = 'off', xlim = c(1, 2))で調整しました。

library(ggdist)
data |> 
  filter(intervention == 3 | intervention == 4) |> 
  filter(occasion == 0) |> 
  mutate(intervention2 = ifelse(intervention == 3,"Gratitude Visit","Recording early memories")) |> 
  ggplot(aes(x = intervention2, y = cesdTotal)) +
  stat_halfeye(width = .5, justification = -.15,point_color = NA, .width = 0) +
  stat_dots(side = "left", dotsize = 1, justification = 1.1, binwidth = 1) +
  geom_boxplot(width = .1,outlier.shape = NA) +
  theme_classic() +
  labs(x = "Intervention", y = "CES-D") +
  coord_cartesian(clip = 'off', xlim = c(1, 2))

• ４群の比較もしてみます。４つなのでifelse()ではなく，case_when()を使っています。

library(ggdist)
data |>  
  filter(occasion == 0) |> 
  mutate(intervention2 = case_when(
    intervention == 1 ~ "Signature Strengths",
    intervention == 2 ~ "Three Good Things",
    intervention == 3 ~ "Gratitude Visit",
    intervention == 4 ~ "Recording early memories")) |> 
  ggplot(aes(x = intervention2, y = cesdTotal)) +
  stat_halfeye(width = .5, justification = -.15,point_color = NA, .width = 0) +
  stat_dots(side = "left", dotsize = 1, justification = 1.1, binwidth = 1) +
  geom_boxplot(width = .1,outlier.shape = NA) +
  theme_classic() +
  labs(x = "Intervention", y = "CES-D") +
  coord_cartesian(clip = 'off', xlim = c(1, 4))

ワーク

• Woodworth et al.(2018)のデータについて，その他の変数に関しても，棒グラフ，ヒストグラム（密度プロット），散布図（バブルプロット），折れ線，箱ひげ図，レインクラウドプロットを描いてみましょう！

２群の平均値差の検定

• データの分布が正規分布：t検定

• サンプルサイズ小さい，データの分布が正規分布から大きく外れる：ウィルコクソン順位和検定

t検定

• データの分布が正規分布の場合t検定を使う。

• Rのt検定では，t.test()を使う。

• まずは平均と標準偏差を並べる。

• t検定は，t.test(従属変数~独立変数, data = data)で実施できるが，以下の引数がある。

  • alternative = c(“two.sided”, “less”, “greater”) : 基本は両側検定

  • mu = 0：もし１標本t検定を行う場合

  • paired = FALSE：対応のあるt検定を行う場合

  • var.equal = FALSE：両群の等分散性が確認できないと通常のt検定はできず， Welch の補正 (Welch’s correlation)が必要になる。rではデフォルトでWelch の補正をかけている。

• まず，必要なデータを整理して，data_2sample_tに保存します。際どいですが，有意差はないようです。ただ，無視できない差に思います。

data_2sample <- data |> 
  filter(intervention == 3 | intervention == 4) |> 
  filter(occasion == 0) |> 
  mutate(intervention = as.factor(intervention)) |> 
  select(intervention, cesdTotal) 

t.test(cesdTotal ~ intervention, 
       data = data_2sample,
       alternative = "two.sided",
       paired = FALSE,
       var.equal = FALSE)

    Welch Two Sample t-test

data:  cesdTotal by intervention
t = 1.8279, df = 139.41, p-value = 0.0697
alternative hypothesis: true difference in means between group 3 and group 4 is not equal to 0
95 percent confidence interval:
 -0.263802  6.727334
sample estimates:
mean in group 3 mean in group 4 
       16.08108        12.84932 

• 効果量も計算します。

cohen.d(data_2sample$cesdTotal~data_2sample$intervention, hedge.correction=TRUE)

Cohen's d

d estimate: 0.3010943 (small)
95 percent confidence interval:
     lower      upper 
-0.0267866  0.6289752 

ウィルコクソン順位和検定

• データに正規性に疑いがある場合，ノンパラメトリック検定のウィルコクソンの順位和検定を用いる。ウィルコクソンの順位和検定では，２群の順位から検定を行う。

• Rのウィルコクソンの順位和検定では，coinパッケージのwilcox_test()を使う。

    Exact Wilcoxon-Mann-Whitney Test

data:  cesdTotal by intervention (3, 4)
Z = 1.4016, p-value = 0.1618
alternative hypothesis: true mu is not equal to 0

２群の比率差の検定

• ２群の比率の検定では，フィッシャーの正確検定を行います。

• まずはクロス集計を計算します。

data_2ratio <- participant |> 
  filter(intervention == 3 | intervention == 4)

table(data_2ratio$sex, data_2ratio$intervention)

   
     3  4
  1 62 61
  2 12 12

• フィッシャーの正確検定fisher.test()で実行できます。

fisher.test(data_2ratio$sex, data_2ratio$intervention)

    Fisher's Exact Test for Count Data

data:  data_2ratio$sex and data_2ratio$intervention
p-value = 1
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.3841867 2.6885100
sample estimates:
odds ratio 
  1.016284 

• カイ二乗検定でも可能です。カイ二乗検定はchisq.test()で実行できます。

chisq.test(data_2ratio$sex, data_2ratio$intervention)

    Pearson's Chi-squared test with Yates' continuity correction

data:  data_2ratio$sex and data_2ratio$intervention
X-squared = 4.1048e-31, df = 1, p-value = 1

３群以上の平均値差の検定

• サンプルサイズが十分に大きい，もしくはデータが正規分布：１要因分散分析

• サンプルサイズが小さい，データが正規分布ではない：クラスカル・ウォリス検定

• もし１要因分散分析やクラスカル・ウォリス検定で有意差があれば，多重比較を行う。

１要因分散分析

• １要因分散分析（一元配置分散分析）は，oneway.test()で実行できる。

• dataのcesdTotalを４群で比較する。まずはデータを用意して，平均と標準偏差を確認する。

data_4group <- data |> 
  filter(occasion == 0) |> 
  mutate(intervention = as.factor(intervention)) |> 
  select(intervention, cesdTotal) 

data_4group |> 
  group_by(intervention) |> 
  summarize(mean = mean(cesdTotal),
            sd = sd(cesdTotal)) |> 
  gt()

• レインクラウドプロットで確認することも重要です。

data |>  
  filter(occasion == 0) |> 
  mutate(intervention2 = case_when(
    intervention == 1 ~ "Signature Strengths",
    intervention == 2 ~ "Three Good Things",
    intervention == 3 ~ "Gratitude Visit",
    intervention == 4 ~ "Recording early memories")) |> 
  ggplot(aes(x = intervention2, y = cesdTotal)) +
  stat_halfeye(width = .5, justification = -.15,point_color = NA, .width = 0) +
  stat_dots(side = "left", dotsize = 1, justification = 1.1, binwidth = 1) +
  geom_boxplot(width = .1,outlier.shape = NA) +
  theme_classic() +
  labs(x = "Intervention", y = "CES-D") +
  coord_cartesian(clip = 'off', xlim = c(1, 4))

• oneway.test()で決定を行います。

oneway.test(cesdTotal ~ intervention, data = data_4group)

    One-way analysis of means (not assuming equal variances)

data:  cesdTotal and intervention
F = 1.8107, num df = 3.00, denom df = 160.69, p-value = 0.1473

• 本来は１要因分散分析で有意差があれば多重比較をしますが，今回は練習のために，有意差はないけど実施します。多重比較はpairwise.t.test()で行います。

pairwise.t.test(data_4group$cesdTotal, data_4group$intervention, p.adjust.method = 'bonferroni')

    Pairwise comparisons using t tests with pooled SD 

data:  data_4group$cesdTotal and data_4group$intervention 

  1    2    3   
2 1.00 -    -   
3 1.00 1.00 -   
4 1.00 0.35 0.42

P value adjustment method: bonferroni 

クラスカル・ウォリス検定

• クラスカル・ウォリス検定は，ウィルコクソンの順位和検定と同様に順位に基づくノンパラメトリック検定です。kruskal_test()で実行できる。

kruskal_test(cesdTotal ~ intervention, data = data_4group)

    Asymptotic Kruskal-Wallis Test

data:  cesdTotal by intervention (1, 2, 3, 4)
chi-squared = 4.6244, df = 3, p-value = 0.2015

• 本来はクラスカル・ウォリス検定で有意差があれば多重比較をしますが，今回は練習のために，有意差はないけど実施します。多重比較はpairwise.wilcox.test()で行います。

pairwise.wilcox.test(data_4group$cesdTotal, data_4group$intervention, p.adjust.method = 'bonferroni')

    Pairwise comparisons using Wilcoxon rank sum test with continuity correction 

data:  data_4group$cesdTotal and data_4group$intervention 

  1    2    3   
2 1.00 -    -   
3 1.00 1.00 -   
4 1.00 0.19 0.97

P value adjustment method: bonferroni 

心理学の文脈で分散分析をする場合

• ANOVA君の利用を勧める。

３群以上の比率差の検定

• 3群以上の比率の差は，カイ二乗検定を用いる。カイ二乗検定は，chisq.test()を用いる。

• dataの介入群ごとの性別の比率をクロス集計で確認する。

table(participant$intervention,participant$sex)

   
     1  2
  1 62 10
  2 66 10
  3 62 12
  4 61 12

• カイ二乗検定を実施する。有意ではない。

chisq.test(participant$intervention,participant$sex)

    Pearson's Chi-squared test

data:  participant$intervention and participant$sex
X-squared = 0.47685, df = 3, p-value = 0.9239

• 本来はカイ二乗検定で有意差があれば多重比較をしますが，今回は練習のために，有意差はないけど実施します。多重比較はpairwise.prop.test()で行います。

pairwise.prop.test(table(participant$intervention,participant$sex), p.adjust.method = 'bonferroni')

    Pairwise comparisons using Pairwise comparison of proportions 

data:  table(participant$intervention, participant$sex) 

  1 2 3
2 1 - -
3 1 1 -
4 1 1 1

P value adjustment method: bonferroni 

相関分析

• ２つの連続変数の関係は散布図や相関係数で表現できる。それが無相関かどうかの帰無仮説検定するのが無相関検定である。

• まずは散布図を確認する。

data |> 
  filter(occasion == 0) |> 
  ggplot(aes(x = cesdTotal, y = ahiTotal)) +
  geom_point() +
  theme_classic() +
  labs(x = "CES-D", y = "AHI")

• 相関係数はcor()で計算できる。相関係数は２つの連続変数間関係についての効果量である。

data_cor <- data |> 
  filter(occasion == 0) 

cor(data_cor$ahiTotal,data_cor$cesdTotal)

[1] -0.7254563

• 無相関検定はcor.test()で実施できる。

cor.test(data_cor$ahiTotal,data_cor$cesdTotal)

    Pearson's product-moment correlation

data:  data_cor$ahiTotal and data_cor$cesdTotal
t = -18.042, df = 293, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.7754145 -0.6664727
sample estimates:
       cor 
-0.7254563 

• 基本的にはcor()やcor.test()で十分だが，多変数の相関関係を確認したい場合は，corrrパッケージが便利かもしれない。

• CES-Dの２０項目の相関関係を検討する。correlate()で相関行列を計算して，cesd_corに格納する。

library(corrr)
cesd_cor <- data_cor |> 
  select(cesd01:cesd20) |> 
  correlate()

• fashion()で見やすくできる。

fashion(cesd_cor)

     term cesd01 cesd02 cesd03 cesd04 cesd05 cesd06 cesd07 cesd08 cesd09 cesd10
1  cesd01           .38    .43   -.21    .41    .39    .37   -.17    .27    .29
2  cesd02    .38           .40   -.23    .20    .35    .29   -.18    .31    .18
3  cesd03    .43    .40          -.44    .43    .73    .59   -.42    .52    .44
4  cesd04   -.21   -.23   -.44          -.26   -.40   -.32    .57   -.33   -.23
5  cesd05    .41    .20    .43   -.26           .45    .46   -.22    .27    .29
6  cesd06    .39    .35    .73   -.40    .45           .67   -.45    .59    .49
7  cesd07    .37    .29    .59   -.32    .46    .67          -.31    .49    .42
8  cesd08   -.17   -.18   -.42    .57   -.22   -.45   -.31          -.40   -.32
9  cesd09    .27    .31    .52   -.33    .27    .59    .49   -.40           .45
10 cesd10    .29    .18    .44   -.23    .29    .49    .42   -.32    .45       
11 cesd11    .28    .25    .34   -.18    .21    .37    .30   -.19    .21    .30
12 cesd12   -.41   -.35   -.64    .49   -.38   -.65   -.53    .56   -.43   -.44
13 cesd13    .33    .28    .40   -.23    .29    .35    .29   -.21    .26    .33
14 cesd14    .26    .31    .45   -.32    .26    .49    .39   -.34    .52    .38
15 cesd15    .23    .09    .16   -.09    .10    .16    .21   -.07    .28    .17
16 cesd16   -.31   -.33   -.60    .48   -.31   -.56   -.43    .53   -.40   -.35
17 cesd17    .27    .29    .49   -.18    .30    .36    .37   -.20    .32    .32
18 cesd18    .47    .33    .69   -.40    .43    .66    .56   -.41    .57    .57
19 cesd19    .24    .17    .35   -.31    .32    .35    .32   -.24    .48    .37
20 cesd20    .35    .32    .52   -.30    .51    .55    .51   -.25    .38    .32
   cesd11 cesd12 cesd13 cesd14 cesd15 cesd16 cesd17 cesd18 cesd19 cesd20
1     .28   -.41    .33    .26    .23   -.31    .27    .47    .24    .35
2     .25   -.35    .28    .31    .09   -.33    .29    .33    .17    .32
3     .34   -.64    .40    .45    .16   -.60    .49    .69    .35    .52
4    -.18    .49   -.23   -.32   -.09    .48   -.18   -.40   -.31   -.30
5     .21   -.38    .29    .26    .10   -.31    .30    .43    .32    .51
6     .37   -.65    .35    .49    .16   -.56    .36    .66    .35    .55
7     .30   -.53    .29    .39    .21   -.43    .37    .56    .32    .51
8    -.19    .56   -.21   -.34   -.07    .53   -.20   -.41   -.24   -.25
9     .21   -.43    .26    .52    .28   -.40    .32    .57    .48    .38
10    .30   -.44    .33    .38    .17   -.35    .32    .57    .37    .32
11          -.33    .19    .21    .16   -.24    .24    .31    .20    .35
12   -.33          -.39   -.46   -.17    .75   -.34   -.63   -.32   -.45
13    .19   -.39           .37    .15   -.32    .28    .40    .26    .33
14    .21   -.46    .37           .26   -.43    .36    .59    .44    .41
15    .16   -.17    .15    .26          -.15    .12    .19    .41    .16
16   -.24    .75   -.32   -.43   -.15          -.29   -.54   -.32   -.43
17    .24   -.34    .28    .36    .12   -.29           .57    .32    .30
18    .31   -.63    .40    .59    .19   -.54    .57           .40    .47
19    .20   -.32    .26    .44    .41   -.32    .32    .40           .36
20    .35   -.45    .33    .41    .16   -.43    .30    .47    .36       

• 見やすそうなプロットもできる。負の相関があるところは逆転項目っぽいので逆転処理してないのかもしれない。

rplot(cesd_cor)

• 相関関係をネットワーク化することもできる。変数間の関係性がわかりやすい。

network_plot(cesd_cor, min_cor = .2)

サンプルサイズ設計

• 論文を書く際にデータ収集前のサンプルサイズ設計が求められるようになっている。

→再現性問題や統計の誤使用の問題からデータ収集中の検定の実施は行うべきではない（有意になるまでNを増やすN増しは絶対にダメ）。データ収集前に必要なサンプルサイズを決めておく必要がある。

→検定力が低い研究は，帰無仮説の棄却ができず結論づけができないので，研究を行う意義が薄くなる。

• 2群比較，相関分析などのサンプルサイズ設計は比較的簡単だが，多変量解析や混合効果モデルなどはかなり難しい。中心的な臨床的疑問をシンプルに絞っておく必要がある。

サンプルサイズの決め方

• 必要なサンプルサイズは，効果量，検定力，有意水準の３つを設定することで計算できる。

• 検定力と有意水準は，慣例に従い，検定力を0.8，有意水準を0.05に設定することが多い。効果量については，研究の対象やアウトカム（効果指標）に応じて設定していく必要がある（文脈によっては小さい効果量に意味があることもある）。

→効果量を決めることが一番むずかしい。効果量が決まれば，サンプルサイズも決まってくる。

効果量の決め方

• Hislop et al. (2014)は，過去の無作為化比較試験で行われたサンプルサイズ設計を系統的に展望して，効果量の設定方法について整理した。

• （１）治療者や患者が実際に差を感じるような意味づけのなされた差を採用する「重要な差(important difference)」（アンカーに基づく方法，分布に基づく方法，医療経済学に基づく方法，標準化効果量に基づく方法）

• （２）過去の研究における群間差を参考にする「現実的な差(realistic difference)」（パイロット研究）

• （３）重要な差と現実的な差の両方を含む決定方法として，意見聴取に基づく方法とエビデンスの展望に基づく方法がある

  →絶対的な基準はないので，研究領域で納得感のある決め方ができれば良い。

pwrで必要サンプルサイズを決める

• pwrパッケージを使えば様々な基本的な分析手法での必要サンプルサイズを計算できる。ここでは，２群のｔ検定と相関分析について実際に計算してみる。

２群の平均値差の必要サンプルサイズ

• 先行研究から２群の平均値差の標準化効果量(Cohenのd)を0.3と見積もったとする。有意水準(0.05), 検出力(0.8)，２群，両側検定とした場合の必要サンプルサイズは以下のように計算する。各群176名以上必要である。

library(pwr)
pwr.t.test(n = NULL, 
           d = 0.3, 
           sig.level = 0.05, 
           power = 0.8, 
           type = "two.sample",
           alternative = "two.sided")

     Two-sample t test power calculation 

              n = 175.3847
              d = 0.3
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

• もし標準化効果量が0.6だったら，以下のように各群45名以上になる。

library(pwr)
pwr.t.test(n = NULL, 
           d = 0.6, 
           sig.level = 0.05, 
           power = 0.8, 
           type = "two.sample",
           alternative = "two.sided")

     Two-sample t test power calculation 

              n = 44.58577
              d = 0.6
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

相関分析の必要サンプルサイズ

• 先行研究から相関係数を0.2と見積もったとする。有意水準(0.05), 検出力(0.8)，両側検定とした場合の必要サンプルサイズは以下のように計算する。各群194名以上必要である。

pwr.r.test(n = NULL, 
           r = 0.2, 
           sig.level = 0.05, 
           power = 0.8,
           alternative = "two.sided")

     approximate correlation power calculation (arctangh transformation) 

              n = 193.0867
              r = 0.2
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

• もし標準化効果量が0.4だったら，以下のように各群46名以上になる。

pwr.r.test(n = NULL, 
           r = 0.4, 
           sig.level = 0.05, 
           power = 0.8,
           alternative = "two.sided")

     approximate correlation power calculation (arctangh transformation) 

              n = 45.91614
              r = 0.4
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

脱落を考慮して調査を計画する

• 必要サンプルサイズは，帰無仮説検定を行う上で必要なサンプルサイズなので，データ収集完了時にこれを下回らないようにする必要がある。

• 脱落を考慮して調査計画を立てれば安心して実施できる。そんなに難しい話ではなく，暫定的に脱落率を決めて（オンライン調査だと30%など），必要サンプルサイズに\(\frac{1}{1-脱落率}\)　をかければ良い。

46*1/(1-0.3)

[1] 65.71429

66名からデータを得るように計画すれば，３割が脱落しても４５名以上からデータが得られる。

ワーク

Woodworth et al.(2018)のデータについて，各種検定を実施してみよう！ご自身の研究についてサンプルサイズ設計をしてみよう！

重回帰

?

回帰でgtsummaryを使うと便利だね。

mod1 <- glm(response ~ trt + age + grade, trial, family = binomial)

t1 <- tbl_regression(mod1, exponentiate = TRUE)

?

?

?

?

?

?

?

?

?

一般化線形モデルによる理解

• 統計解析は種類が多い…

一般化線形モデルによる理解

• これらの解析は，線形モデル（一般線形，一般化線形）として理解できる。

一般化線形モデルによる理解

一般化線形モデルによる理解

一般化線形モデルによる解析手順

1.モデルの特定

→結果変数と説明変数を連結する関係式＆反応変数の確率分布の２つを指定

2.モデルのパラメータの推定

3.モデルの妥当性のチェック

→モデルはどのくらいよくデータに当てはまり、データを要約しているか

4.推測

→モデルのパラメータに関する仮説検定と信頼区間の計算および結果の解釈

一般線形モデルによる解析手順

1. モデルの特定：結果変数と説明変数を連結する関係式＆反応変数の確率分布の２つを指定

• 一般線形モデルにおける前提は以下の４つ

一般線形モデルにおける前提：正規性、等分散性、線形性

一般線形モデルにおける前提：独立性(個々の観測値が独立)

• 以下のようなデータを抽出した場所や反復測定する個人で階層構造をもつデータは独立性のないデータになる。

→独立性のないデータは、一般化線形混合モデルを用いる。

一般線形モデルによる分散分析の理解

• t検定、分散分析、重回帰は、一般線形モデルから全て以下のような回帰モデルを作ることができる。

正規性・等分散性が満たされない場合

• ヒストグラム、歪度、尖度、Kolmogorov-Smirnov（コロモゴロフ・スミノフ）検定などで確認

• データの分布によっては対数変換、外れ値の吟味

→現実にあわないデータ処理によって、データが歪む可能性もあるので注意

• 平均値差の検定の場合は、Welch補正を使う。一般線形モデルでは、ロバスト回帰分析もある。

• 一般線形モデルの仮定が難しいデータの場合は、ノンパラメトリック分析へ

連続変数のパラメトリック解析：一般(化)線形混合モデル

一般化線形混合モデル

線形モデルからランダム切片モデルへ

ランダム傾きモデルとランダム切片＆傾きモデル

• ランダム傾きモデル：切片は全体と同じだが、傾きが個人によって異なる

• ランダム切片＆傾きモデル：切片も傾きも個人によって異なる。

一般化線形混合モデルのメリット

• これまで使用されてきた反復測定分散分析に比べて・・・

1. ランダムな欠測には影響をうけない（利用可能なデータをすべて使える）

2. 不等分散や非独立性に強い

3. 少サンプルでもバイアスが少ない（正確さが増すので検出力も上がる）

4. 査読で求められることが増えてきている・・・

?

はしらないようにしてある。

library(lme4)
library(phia)
library(lmerTest)
BtheBlong %>% 
  dplyr::group_by(treatment,time) %>% 
  dplyr::summarise(Mean=mean(bdi, na.rm=T), SD=sd(bdi, na.rm=T), Min=min(bdi, na.rm=T), Max=max(bdi, na.rm=T))


BtheBlong %>% 
  ggplot(aes(x=time, y=bdi, colour=treatment, group=id))+
  geom_line(alpha=.5)+
  scale_x_discrete(label=c("Post","1mFU","3mFU","6mFU"),name="Time")+
  scale_y_continuous(name="BDI")


#ランダム切片モデルを使った一般化線形混合モデル解析を実施した（実際は複数のモデルで比較する必要がある）。
result <- lmer(bdi~bdi.pre + time*treatment+(1|id),data = BtheBlong, REML = FALSE,na.action = na.omit)
summary(result)

result95ci <- confint(result,method = "Wald")
result95ci

resultAdjM <- interactionMeans(result)
print(resultAdjM)

resultAdjM %>% 
  dplyr::rename(adjMean = `adjusted mean`, se = `SE of link`) %>% 
  ggplot(aes(x = time, y = adjMean, colour = treatment, group = treatment)) +
  geom_line(alpha=.5)+theme_classic(base_size = 12,base_family = "") +
  geom_errorbar(aes(ymin = adjMean - 1.96*se,　ymax = adjMean + 1.96*se,width = 0.1))+
  scale_x_discrete(label=c("Post","1mFU","3mFU","6mFU"),name="Time")+
  scale_y_continuous(name="Adjusted Mean") +
  labs(colour = "Treatment")

?

?

ワーク

?

構造方程式モデル

?

媒介分析

?

確証的因子分析

?

?

?

?

?

?

?

?

?

?

?

?

奥村さんのコード

#----------------------------#
#á@ÉfÅ[É^ÇÃì«Ç›çûÇ›
#----------------------------#
##(1) égópÉâÉCÉuÉâÉäÇÃèÄîı
install.packages(c("psych", "paran"), contriburl=contrib.url("http://cran.md.tsukuba.ac.jp/"))  #1âÒÉCÉìÉXÉgÅ[ÉãÇ∑ÇÍÇŒ2ìxñ⁄ÇÕïsóv
library(psych)
library(paran)

##(2) ì«Ç›çûÇ›
setwd("W:/workshop2011/archives/04/")
dat <- read.table("HolzingerGW.csv", sep=",", header=T)
head(dat)
tail(dat)


#----------------------------#
#áAàÍïœó âêÕÇ∆ìÒïœó âêÕ
#----------------------------#
summary(dat)
describe(dat)
cor(dat)

boxplot(dat)
graphics.off()


#----------------------------#
#áBàˆéqêîÇÃåàíË
#----------------------------#

##(1) ÉJÉCÉUÅ[äÓèÄÇ∆ÉXÉNÉäÅ[äÓèÄ
x <- eigen(cor(dat))
plot(x$values)
abline(h=1)
graphics.off()


##(2) ïΩçsï™êÕ
paran(dat, cfa=TRUE, all=TRUE, graph=TRUE, iterations=1000)
graphics.off()


##(3) ç≈è¨ïΩãœïŒëää÷
VSS(dat, n= 4, plot=FALSE)


##(4) BICÇ∆RMSEA
fit1 <- fa(r=dat, nfactors=1, fm="ml", rotate="promax")
fit2 <- fa(r=dat, nfactors=2, fm="ml", rotate="promax")
fit3 <- fa(r=dat, nfactors=3, fm="ml", rotate="promax")

rbind(fit1$RMSEA, fit2$RMSEA, fit3$RMSEA)
rbind(fit1$BIC, fit2$BIC, fit3$BIC)


##(5) ã§í ê´ÅCàˆéqç\ë¢ÅCàˆéqÉpÉ^ÉìÅCàˆéqä‘ëää÷
fit2
fit2$communality                  #ã§í ê´
fit2$Structure                    #àˆéqç\ë¢
print(fit2$loadings, cutoff=0)    #àˆéqÉpÉ^Éì
fit2$Phi                          #àˆéqä‘ëää÷

?

?

メタ分析の一般的手順

系統的展望のポイントについては，１日目の批判的吟味を参照ください。ここでは，解析部分のメタ分析に焦点を当てます。

1. 個々の研究の効果量を計算

2. 個々の研究の効果量を統合

3. 統計学的異質性を評価

4. 事前に計画したサブグループ解析

5. 感度分析

6. 非報告バイアス（出版バイアス）の検討

連続変数における２群の差の効果量

２群の差は，標準化平均値差(Hedgesのg)を用いる

• \(Hedgesのg = \frac{介入群の平均-統制群の平均}{２群をプールした標準偏差}\)

• \(2群をプールした標準偏差=\sqrt{\frac{(n_{介入}-1)\hat{\sigma}^{2}_{介入}+(n_{統制}-1)\hat{\sigma}^{2}_{統制}}{(n_{介入}+n_{統制}-2)}}\)

\(\hat{\sigma}^{2}\)は不偏分散。他にCohenのdやGlassの⊿もあるが，バイアスを修正した効果量のHedgesのgがよく使用される

• \(標準化平均値差の分散=\frac{n_{介入}+n_{統制}}{n_{介入}n_{統制}}+\frac{g^2}{2(n_{介入}+n_{統制})}\)

• \(標準化平均値差の標準誤差=\sqrt{標準化平均値差の分散}\)

分散と標準偏差は効果量とサンプルサイズで決まる

２つの連続変数の関連の効果量

２変数間の関連の強さは，ピアソンの積率相関係数を用いる

• \(r = \frac{共分散_{x,y}}{標準偏差_{x}標準偏差_{y}}\)

• \(相関係数の分散 = \frac{(1-r^2)^2}{n-1}\),\(相関係数のSE = \frac{1-r^2}{\sqrt{n-1}}\)

• 統計学的特性の望ましさから(サンプルサイズが小さいとrのseにはバイアスが生じる)，フィッシャーのZ変換を用いることが多い

• \(z = 0.5ln(\frac{1+r}{1-r})\), \(zの分散 = \frac{1}{n-3}\), \(zのSE = \frac{1}{\sqrt{n-3}}\)

• zからrへは\(\frac{exp(2z)-1}{exp(2z)+1}\)で変換できる(信頼区間も同様)。

２値変数（あり・なし）の効果量

1. リスク差=介入群のありの比率-統制群のありの比率

2. \(治療必要数=\frac{1}{リスク差}\)

3. \(リスク比=\frac{介入群のありの比率}{統制群のありの比率}\)

4. \(オッズ比=\frac{\frac{介入群のありの比率}{1-介入群のありの比率}}{\frac{統制群のありの比率}{1-統制群のありの比率}}\)

リスク差や治療必要数(Number needed to treat)は解釈しやすいが，計算上オッズ比が便利なので，オッズが使用される。さらに対数オッズ比(ln(オッズ比))を使う（1を起点に解釈→0を起点に解釈，分散も計算がしやすい）。

オッズ比の分散=介入群あり比率+介入群なし比率+統制群あり比率+統制群なし比率

固定効果モデル

\(統合された効果量=\frac{\Sigma W_{i}Y_{i}}{\Sigma W_{i}}=\frac{(各研究の重みW_{i}*各研究の効果量Y_{i})の合計}{各研究の重みW_{i}の合計}\)

\(W_{i} = \frac{1}{V_{i}} = \frac{1}{効果量の分散V_{i}}\)

固定効果モデル

• 統合された効果量の標準誤差は，各研究の重みの合計から計算される

\(統合された効果量の標準誤差=\sqrt{\frac{1}{各研究の重みW_{i}の合計}}\)

• 95%信頼区間(CI)は，以下のように計算する。

\(95\%CIの下限=統合された効果量-1.96*統合された効果量の標準誤差\)

\(95\%CIの上限=統合された効果量+1.96*統合された効果量の標準誤差\)

ランダム効果モデル

ランダム効果モデル

• 平均値差の統合(DerSimonian-Laird法)では，以下のように，重み付けに研究間分散(\(\tau^2\))と研究内分散(\(V_{i}\)）を用いる。

1.固定効果モデルを用いて統合された効果量を算出し，Qを計算する(Tは個々の研究の効果量\(\bar{T}\)は統合された効果量)。Qは研究間のばらつきの程度。

\(Q = \Sigma W_{i}(T_{i}-\bar{T})^2\)

2.Qを用いて，研究間分散(\(\tau^2\))を計算する(kはメタ分析に含めた研究数)

\(\tau^2 = \frac{Q-(k-1)}{\Sigma W_{i} - \frac{\Sigma W^2_{i}}{\Sigma W_{i}}}\)

3.重み付けの計算に\(\tau^2\)\(V_{i}\)を用いる以外は，以降は固定効果モデルと同じ。

\(W_{i} = \frac{1}{V_{i}+\tau^2}\)

ランダム効果モデル

• DerSimonian-Lairdの他に，Paule-Mandel, Restricted Maximum-Likelihood, Maximum-likelihood, Hunter-Schmidt, Sidik-Jonkman, Hedges, Empirical Bayesがある(Rのmeta パッケージで利用可能)。

• DerSimonian-Lairdがよく使われるが，それはRevManで使われているなどの理由であって，Maximum-Likelihood, Sidik-Jonkman, Empirical Bayesの方が研究間変動の推定には良い性質がある(Harrer rt al., 2019)

• ランダム効果か固定効果かは事前登録の段階で決める。

概念的異質性

• 論文を精査することで明らかとなる質的な研究間の違い

• 臨床的異質性：年齢，性別，人種，重症度，合併症など患者集団の違い（心理学だと対象者属性の異質性）

• 方法論的異質性：治療法，治療期間，治療者の経験など方法の違い（心理学だと測定法，実験操作などの違い）

統計学的異質性

• 統計学的に明らかとなる研究ごとの効果のばらつき

→統計学的異質性も概念的異質性もある場合：(1)メタ分析しない,(2)サブグループ解析・メタ回帰分析の実施

統計学的異質性

• ランダム効果モデルのQを使った検定もあるが，研究数が多いと有意になりやすい（有意じゃなくても検定力低いだけかも）

→ 検定ではなく，Higgins et al.(2003)の\(I^2=\frac{Q-(k-1)}{Q}100\)のような記述的指標を用いて判断をする。 \(I^2\)は，効果量のバラツキの内，効果量の異質性が占める割合を％で表している。

• \(I^2\)のおおまかな目安

  • 40%以下 : 異質性は問題ではないかも

  • 30~60% : 中程度の異質性があるかも

  • 50~90% : 実際に異質性があるかも

  • 75~100% : 無視できない異質性がある

サブグループ解析

• 異質性を生むと考えられるサブグループを特定して（患者の特性，方法など），分けて効果量の統合をする。

→プロトコルの段階で計画を立てる

感度分析

• 適格基準や統合方法を変えても今回の結果が頑健に得られるか確認する。

→(1)適格基準の変更（患者の特徴を変えても変わらない？）,(2)分析データの変更（調整済み効果量では？），(3)分析法の変更(固定効果とランダム効果モデルで大きな差はない？)

出版バイアス

MASEM(Meta-analytic structural equation modeling)

• メタ分析をSEMの枠組みで実行することができる。

• SEMを活用することで，多変量メタ分析ができる。例えば，臨床試験には複数のアウトカムが設定されているが，アウトカム間に関連があることもあり，それを考慮した上で効果を検討することができる。

• メタ分析によって複数の研究で報告されている相関行列を統合し，SEMで確証的因子分析を実施したり，媒介分析を実施したりもでる。

• 詳しくは，Jak(2015)，Cheung(2015)，宇佐美(2018)を参照

• RのmetaSEMパッケージが使える。

Rでメタ分析

• Rのメタ分析の資料としては，『Doing Meta-Analysis in R : A Hands-on Guide』(Harrer, M., Cuijpers, P., Furukawa,T.A., Ebert, D.D.)を超える資料はないのではないかと思います。無料で読めます。

2群の平均値差:データ準備

# パッケージの読みこみ
library(dmetar)
library(meta)
# データを作る（Mは平均，Sは標準偏差，Nは人数，eが介入・実験群，cが統制群）
Me <- c(35, 40, 35, 22, 38, 31, 29, 28)
Se <- c(7.5, 8.2, 5.9, 10.4, 8.8, 9.1, 6.4, 4.6)
Ne <- c(15, 32, 18, 22, 27, 42, 29, 17)
Mc <- c(20, 23, 30, 10, 22, 25, 19, 20)
Sc <- c(6.5, 4.3, 5.7, 7.8, 6.9, 10.4, 7.1, 6.5)
Nc <- c(16, 30, 19, 21, 29, 40, 32, 18)
Author<- c("Senshu (1923)","Waseda (1920)","Doshisha (1920)","Kangaku (1932)","Keio (1920)","Ritsumei (1922)","Hyokyo (1978)","Toyo (1928)")
Subgroup <- c("East","East","West","West","East","West","West","East")
# データフレームを作る
data_meta <- data.frame(Author, Me, Se, Ne, Mc, Sc, Nc, Subgroup)

2群の平均値差：メタ分析の実施

• metaパッケージの連続変数用メタ分析関数のmetacontを使う。

• comb.fixedが固定効果，comb.randomがランダム効果になる。

• method.smdは，デフォルトが”Hedges”だが，“Cohen”や”Glass”も選べる。

• method.tauは，デフォルトがDL(DerSimonian-Laird)だが他にも選択できる。今回は，最尤法(ML)を選択した。

result_meta <- metacont(Ne,Me,Se,Nc,Mc,Sc,
                  data = data_meta,
                  studlab = paste(Author),
                  comb.fixed = TRUE,
                  comb.random = TRUE,
                  prediction = FALSE,
                  sm = "SMD",
                  method.smd = "Hedges",
                  method.tau = "ML")

2群の平均値差：メタ分析結果

result_meta

Number of studies: k = 8
Number of observations: o = 407

                        SMD           95%-CI     z  p-value
Common effect model  1.3693 [1.1470; 1.5917] 12.07 < 0.0001
Random effects model 1.4897 [1.0553; 1.9240]  6.72 < 0.0001

Quantifying heterogeneity:
 tau^2 = 0.2783 [0.0839; 1.6462]; tau = 0.5276 [0.2897; 1.2831]
 I^2 = 77.7% [55.8%; 88.7%]; H = 2.12 [1.50; 2.97]

Test of heterogeneity:
     Q d.f.  p-value
 31.33    7 < 0.0001

Details on meta-analytical method:
- Inverse variance method
- Maximum-likelihood estimator for tau^2
- Q-Profile method for confidence interval of tau^2 and tau
- Hedges' g (bias corrected standardised mean difference; using exact formulae)

2群の平均値差：フォレストプロット

forest(result_meta)

2群の平均値差：異質性の確認

summary(result_meta)

                   SMD           95%-CI %W(common) %W(random)
Senshu (1923)   2.0867 [1.1906; 2.9828]        6.2       10.1
Waseda (1920)   2.5400 [1.8611; 3.2188]       10.7       12.3
Doshisha (1920) 0.8437 [0.1679; 1.5195]       10.8       12.4
Kangaku (1932)  1.2770 [0.6157; 1.9383]       11.3       12.5
Keio (1920)     2.0041 [1.3542; 2.6540]       11.7       12.6
Ritsumei (1922) 0.6093 [0.1658; 1.0527]       25.1       14.9
Hyokyo (1978)   1.4568 [0.8878; 2.0257]       15.3       13.5
Toyo (1928)     1.3813 [0.6352; 2.1273]        8.9       11.6

Number of studies: k = 8
Number of observations: o = 407

                        SMD           95%-CI     z  p-value
Common effect model  1.3693 [1.1470; 1.5917] 12.07 < 0.0001
Random effects model 1.4897 [1.0553; 1.9240]  6.72 < 0.0001

Quantifying heterogeneity:
 tau^2 = 0.2783 [0.0839; 1.6462]; tau = 0.5276 [0.2897; 1.2831]
 I^2 = 77.7% [55.8%; 88.7%]; H = 2.12 [1.50; 2.97]

Test of heterogeneity:
     Q d.f.  p-value
 31.33    7 < 0.0001

Details on meta-analytical method:
- Inverse variance method
- Maximum-likelihood estimator for tau^2
- Q-Profile method for confidence interval of tau^2 and tau
- Hedges' g (bias corrected standardised mean difference; using exact formulae)

異質性に関連した外れ値の検討，Influence Analyses，GOSH Plot Analysisは省略しました。詳細は，Doing Meta-Analysis in R : A Hands-on Guideを確認ください

2群の平均値差：出版バイアス

funnel(result_meta,　xlab = "g",studlab = TRUE)

2群の平均値差：出版バイアス

dmetarパッケージのeggers.testを使う。

eggers.test(x = result_meta)

Eggers' test of the intercept 
============================= 

 intercept       95% CI     t          p
     6.568 0.48 - 12.65 2.116 0.07877153

Eggers' test does not indicate the presence of funnel plot asymmetry. 

2群の平均値差：出版バイアスの調整(trim and fill法)

trimfill(result_meta)

Number of studies: k = 10 (with 2 added studies)
Number of observations: o = 500

                        SMD           95%-CI    z  p-value
Random effects model 1.2193 [0.7183; 1.7204] 4.77 < 0.0001

Quantifying heterogeneity:
 tau^2 = 0.5286 [0.2202; 2.3463]; tau = 0.7271 [0.4692; 1.5318]
 I^2 = 83.5% [71.1%; 90.5%]; H = 2.46 [1.86; 3.25]

Test of heterogeneity:
     Q d.f.  p-value
 54.39    9 < 0.0001

Details on meta-analytical method:
- Inverse variance method
- Maximum-likelihood estimator for tau^2
- Q-Profile method for confidence interval of tau^2 and tau
- Trim-and-fill method to adjust for funnel plot asymmetry (L-estimator)

2群の平均値差：出版バイアスの調整結果

funnel(trimfill(result_meta),xlab = "g",studlab = TRUE)

今回触れていないこと

• サンプルサイズを考慮して真の効果量を検討するP-curve分析(Simonsohn et al., 2014)

• Risk of　Biasのプロット方法

• サブグループ解析

• ２値変数，相関係数のメタ分析

• ベイジアンメタ分析

• ネットワークメタ分析

• メルチレベルメタ分析

Doing Meta-Analysis in R : A Hands-on Guideでより詳しく学ぶことができます。

レポート２

関心のあるテーマでオープンデータ・模擬データを探して，Quartoを使って解析レポートを作成しましょう。